#include "freeboundary.h"

struct particle_line *init_particle(	char *problem,
					int *NP,
					int ppc,
					double dx,
					double dy,
					int imax,
					int jmax,
					double **U,
					double **V,
					int *npart)
{
	double height, width;
	double radius;
	double mp_x, mp_y;
	double x_par, y_par;
	double x_pos, y_pos;
	double x_right, x_left;
	double y_top, y_bottom;
	double y_start;
	double rad_square;

	struct particle_line *list;
	struct particle *pnew;
	struct particle *prev;

	x_par=dx/ppc;
	y_par=dy/ppc;

	if (strcmp(problem, "Breaking_Dam")==0)
	{
		*NP = 1;
		list = (struct particle_line *)malloc((*NP)*(sizeof(struct particle_line)));
		width = imax*(dx/5);
		list[0].length = 0;

		x_pos = x_par/2;
		while(x_pos<=width)
		{
			y_pos = y_par/2;

			while(y_pos<=(jmax*dy))
   			{
				pnew = (struct particle *)malloc(sizeof(struct particle));

				if(list[0].length==0)
				{
					list[0].Particle = pnew;
				}
				else
				{
				 	prev->next=pnew;
        		}

				list[0].length++;

				pnew->x = x_pos;
				pnew->y = y_pos;

				pnew->next = NULL;
				prev = pnew;

				y_pos = y_pos + y_par;
   			}
			x_pos = x_pos + x_par;
		}
		*npart = list[0].length;
	}


	if (strcmp(problem, "Falling_Drops")==0)
	{
		*NP = 2;
		list = (struct particle_line *)malloc((*NP)*(sizeof(struct particle_line)));
		height=jmax*dy/3;
		radius=jmax*dy/10;
		mp_x=imax*dx/2;
		mp_y=(2.0/3.0)*(jmax*dy);

		list[0].length = 0;
		list[1].length = 0;

		x_pos = x_par/2;

		while(x_pos<=(imax*dx))
		{
			y_pos = y_par/2;

			while(y_pos<=height)
   			{
				pnew = (struct particle *)malloc(sizeof(struct particle));

				if(list[0].length==0)
				{
					list[0].Particle = pnew;
				}
				else
				{
				 	prev->next=pnew;
        		}

				list[0].length++;

				pnew->x = x_pos;
				pnew->y = y_pos;

				pnew->next = NULL;
				prev = pnew;

				y_pos = y_pos + y_par;
   			}
			x_pos = x_pos + x_par;
		}

		x_left = mp_x - radius;
		y_bottom = mp_y - radius;
		x_right = mp_x + radius;
		y_top = mp_y + radius;

		x_pos = x_par/2;
		y_pos = y_par/2;


		while(x_pos<=x_left) {
			x_pos = x_pos + x_par;
		}

        	while(y_pos<=y_bottom) {
        		y_pos = y_pos + y_par;
		}

		y_start = y_pos;


		while(x_pos<=x_right) {
			y_pos = y_start;

			while(y_pos<=y_top) {
				rad_square = pow((x_pos-mp_x),2)+pow((y_pos-mp_y),2);

	        		if(rad_square <= pow(radius,2)) {
          				pnew = (struct particle *)malloc(sizeof(struct particle));

					if(list[1].length == 0) {
						list[1].Particle = pnew;
					}
					else {
				 		prev->next=pnew;
        				}

					list[1].length++;
					pnew->x = x_pos;
					pnew->y = y_pos;

					pnew->next = NULL;
					prev = pnew;
				}
				y_pos = y_pos + y_par;
			}
			x_pos = x_pos + x_par;
		}
		*npart = list[0].length + list[1].length;

	}
	return list;
}

void mark_cells(int **Flag,
		double dx,
		double dy,
		int imax,
		int jmax,
		int N,
		struct particle_line *Particles_Lines)
{
	int i,j,k;
	double x_pos,y_pos;
	struct particle *ptr;

	init_imatrix( Flag, 0, imax + 1, 0 , jmax + 1, 0);

	/* Set the Fluid cells */
	for ( k = 0; k < N; k++ ) {
		ptr = Particles_Lines[k].Particle;

		while(ptr!=NULL) {
			x_pos = ptr -> x;
			y_pos = ptr -> y;

			i = (int)(x_pos/dx) + 1;
			j = (int)(y_pos/dy) + 1;

			Flag[i][j] = C_F;
			ptr = ptr -> next;
		}
	}

	for ( i = 1; i <= imax; i++ ) {
		for ( j = 1; j <= jmax; j++ ) {
			if (!is_fluid(Flag[i][j]))
				Flag[i][j] = C_Y;
		}
	}

	for ( i = 1; i <= imax; i++ ) {
		for ( j = 1; j <= jmax; j++ ) {
			if (is_fluid(Flag[i][j])) {
	 			if (is_empty(Flag[i][j+1]))
	 				Flag[i][j]|=Y_N;
				if (is_empty(Flag[i][j-1]))
					Flag[i][j]|=Y_S;
				if (is_empty(Flag[i-1][j]))
					Flag[i][j]|=Y_W;
				if (is_empty(Flag[i+1][j]))
					Flag[i][j]|=Y_E;
			}
		}
	}
}


void set_uvp_rand (
			double dt,
			double dx,
			double dy,
			int imax,
			int jmax,
			double Re,
			double **U,
			double **V,
			double **P,
			int **Flag,
			double GX,
			double GY
		)
{
	int i, j;

	/* 1st pass: calculate U & V */
	for ( i = 1; i <= imax; i++ ) {
		for ( j = 1; j <= jmax; j++ ) {
			switch ( Flag[i][j] ) {
				/* Cells with 1 empty neighbour */
				case C_E:
					U[i][j] = U[i - 1][j] - (dx/dy)*(V[i][j] - V[i][j - 1]);
					if (is_empty(Flag[i + 1][j - 1])) {
						V[i + 1][j - 1] = V[i][j - 1] - \
								(dx/dy)*(U[i][j] - U[i][j - 1]);
					}
					break;
				case C_N:
					/* TODO: Check the 1st term of U */
					V[i][j] = V[i][j - 1] - (dy/dx)*(U[i][j] - U[i - 1][j]);
					if (is_empty(Flag[i - 1][j + 1])) {
						U[i - 1][j + 1] = U[i - 1][j] - \
								(dy/dx)*(V[i][j] - V[i - 1][j]);
					}
					break;
				case C_W:
					U[i - 1][j] = U[i][j] + (dx/dy)*(V[i][j] - V[i][j - 1]);
					if (is_empty(Flag[i - 1][j - 1])) {
						V[i - 1][j - 1] = V[i][j - 1] + \
								(dx/dy)*(U[i - 1][j] - U[i - 1][j -1]);
					}
					break;
				case C_S:
					V[i][j - 1] = V[i][j] + (dy/dx)*(U[i][j] - U[i - 1][j]);
					if (is_empty(Flag[i - 1][j - 1])) {
						U[i - 1][j - 1] = U[i - 1][j] + \
								(dy/dx)*(V[i][j - 1] - V[i - 1][j - 1]);
					}
					break;
				/* Cells with 2 empty neighbours sharing a common corner */
				case C_NE:
					U[i][j] = U[i - 1][j];
					V[i][j] = V[i][j - 1];
					if (is_empty(Flag[i + 1][j - 1])) {
						V[i + 1][j - 1] = V[i][j - 1] - \
								(dx/dy)*(U[i][j] - U[i][j - 1]);
					}
					if (is_empty(Flag[i + 1][j + 1])) {
						U[i][j + 1] = U[i][j];
						V[i + 1][j] = V[i][j];
					}
					if (is_empty(Flag[i - 1][j + 1])) {
						U[i - 1][j + 1]	= U[i - 1][j] - \
								(dy/dx)*(V[i][j] - V[i - 1][j]);
					}

					break;
				case C_SW:
					U[i - 1][j] = U[i][j];
					V[i][j - 1] = V[i][j];
					if (is_empty(Flag[i - 1][j - 1])) {
						U[i - 1][j - 1] = U[i - 1][j];
						V[i - 1][j - 1] = V[i][j - 1];
					}
					break;
				case C_NW:
					U[i - 1][j] = U[i][j];
					V[i][j] = V[i][j - 1];
					if (is_empty(Flag[i - 1][j + 1])) {
						U[i - 1][j + 1] = U[i - 1][j];
						V[i - 1][j] = V[i][j];
					}
					if (is_empty(Flag[i - 1][j - 1])) {
						V[i - 1][j - 1] = V[i][j - 1] + \
							 	(dx/dy)*(U[i - 1][j] - U[i - 1][j - 1]);
					}
					break;
				case C_SE:
					U[i][j] = U[i - 1][j];
					V[i][j - 1] = V[i][j];

					if (is_empty(Flag[i + 1][j - 1])) {
						U[i][j - 1] = U[i][j];
						V[i + 1][j - 1] = V[i][j - 1];
					}
					if (is_empty(Flag[i - 1][j - 1])) {
						U[i - 1][j - 1] = U[i - 1][j] + \
								(dy/dx)*(V[i][j - 1] - V[i - 1][j - 1]);
					}
					break;
				/* Cells with 2 opposite empty neighbours */
				case C_WE:
					U[i][j] = U[i][j] + dt*GX;
					U[i - 1][j] = U[i - 1][j] + dt*GX;

					if (is_empty(Flag[i - 1][j - 1])) {
						V[i - 1][j - 1] = V[i][j - 1] + \
								(dx/dy)*(U[i - 1][j] - U[i - 1][j - 1]);
					}
					if (is_empty(Flag[i + 1][j - 1])) {
						V[i + 1][j - 1] = V[i][j - 1] - \
								(dx/dy)*(U[i][j] - U[i][j - 1]);
					}
					break;
				case C_NS:
					V[i][j] = V[i][j] + dt*GY;
					V[i][j - 1] = V[i][j - 1] + dt*GY;

					if (is_empty(Flag[i - 1][j + 1])) {
						U[i - 1][j + 1] = U[i - 1][j] - \
								(dy/dx)*(V[i][j] - V[i - 1][j]);
					}
					if (is_empty(Flag[i - 1][j - 1])) {
						U[i - 1][j - 1] = U[i - 1][j] + \
								(dy/dx)*(V[i][j - 1] - V[i - 1][j - 1]);
					}
					break;
				/* Cells with 3 empty neighbours */
				case C_NWE:
					U[i][j] = U[i][j] + dt*GX;
					U[i - 1][j] = U[i - 1][j] + dt*GX;
					V[i][j] = V[i][j - 1] - (dy/dx)*(U[i][j] - U[i - 1][j]);

					if (is_empty(Flag[i - 1][j + 1])) {
						U[i - 1][j + 1] = U[i - 1][j];
						V[i - 1][j] = V[i][j];
					}
					if (is_empty(Flag[i + 1][j + 1])) {
						U[i][j + 1] = U[i][j];
						V[i + 1][j] = V[i][j];
					}
					if (is_empty(Flag[i - 1][j - 1])) {
						V[i - 1][j - 1] = V[i][j - 1] + \
								(dx/dy)*(U[i - 1][j] - U[i - 1][j - 1]);
					}
					if (is_empty(Flag[i + 1][j - 1])) {
						V[i + 1][j - 1] = V[i][j - 1] - \
								(dx/dy)*(U[i][j] - U[i][j - 1]);
					}
					break;
				case C_SWE:
					U[i][j] = U[i][j] + dt*GX;
					U[i - 1][j] = U[i - 1][j] + dt*GX;
					V[i][j - 1] = V[i][j] + (dy/dx)*(U[i][j] - U[i - 1][j]);

					if (is_empty(Flag[i - 1][j - 1])) {
						U[i - 1][j - 1] = U[i - 1][j];
						V[i - 1][j - 1] = V[i][j - 1];
					}
					if (is_empty(Flag[i + 1][j - 1])) {
						U[i][j - 1] = U[i][j];
						V[i + 1][j - 1] = V[i][j - 1];
					}
					break;
				case C_NSW:
					V[i][j] = V[i][j] + dt*GY;
					V[i][j - 1] = V[i][j - 1] + dt*GY;
					U[i - 1][j] = U[i][j] + (dx/dy)*(V[i][j] - V[i][j - 1]);

					if (is_empty(Flag[i - 1][j - 1])) {
						U[i - 1][j - 1] = U[i - 1][j];
						V[i - 1][j - 1] = V[i][j - 1];
					}
					if (is_empty(Flag[i - 1][j + 1])) {
						V[i - 1][j] = V[i][j];
						U[i - 1][j + 1] = U[i - 1][j];
					}
					break;
				case C_NSE:
					V[i][j] = V[i][j] + dt*GY;
					V[i][j - 1] = V[i][j - 1] + dt*GY;
					U[i][j] = U[i - 1][j] - (dx/dy)*(V[i][j] - V[i][j - 1]);

					if (is_empty(Flag[i + 1][j - 1])) {
						U[i][j - 1] = U[i][j];
						V[i + 1][j - 1] = V[i][j - 1];
					}

					if (is_empty(Flag[i + 1][j + 1])) {
						U[i][j + 1] = U[i][j];
						V[i + 1][j] = V[i][j];
					}
					if (is_empty(Flag[i - 1][j - 1])) {
						U[i - 1][j - 1] = U[i - 1][j] + \
								(dy/dx)*(V[i][j - 1] - V[i - 1][j - 1]);
					}
					if (is_empty(Flag[i - 1][j + 1])) {
						U[i - 1][j + 1] = U[i - 1][j] - \
								(dy/dx)*(V[i][j] - V[i - 1][j]);
					}
					break;
				/* Isolated cell with 4 neighbours */
				case C_NSWE:
					V[i][j] = V[i][j] + dt*GY;
					V[i][j - 1] = V[i][j - 1] + dt*GY;
					U[i][j] = U[i][j] + dt*GX;
					U[i - 1][j] = U[i - 1][j] + dt*GX;

					if (is_empty(Flag[i + 1][j - 1])){
						U[i][j - 1] = U[i][j];
						V[i + 1][j - 1] = V[i][j - 1];
					}
					if (is_empty(Flag[i + 1][j + 1])) {
						U[i][j + 1] = U[i][j];
						V[i + 1][j] = V[i][j];
					}
					if (is_empty(Flag[i - 1][j - 1])) {
						U[i - 1][j - 1] = U[i - 1][j];
						V[i - 1][j - 1] = V[i][j - 1];
					}
					if (is_empty(Flag[i - 1][j + 1])) {
						U[i - 1][j + 1] = U[i - 1][j];
						V[i - 1][j] = V[i][j];
					}
					break;
			}
		}
	}

	/* 2nd pass: calculate P */
	for ( i = 1; i <= imax; i++ ) {
		for ( j = 1; j <= jmax; j++ ) {
			switch (Flag[i][j]) {
				/* Cells with 1 empty neighbour */
				case C_E:
					P[i][j] = (2.0/Re)*(U[i][j] - U[i - 1][j])/dx;
					break;
				case C_N:
					P[i][j] = (2.0/Re)*(V[i][j] - V[i][j - 1])/dy;
					break;
				case C_W:
					P[i][j] = (2.0/Re)*(U[i][j] - U[i - 1][j])/dx;
					break;
				case C_S:
					P[i][j] = (2.0/Re)*(V[i][j] - V[i][j - 1])/dy;
					break;
				/* Cells with 2 neighbours sharing a common corner */
				case C_NE:
					P[i][j] = (1.0/(2.0*Re))*((U[i][j] + U[i - 1][j] - U[i][j - 1] - U[i - 1][j - 1])/dy + \
						(V[i][j] + V[i][j - 1] - V[i - 1][j] - V[i - 1][j - 1])/dx);
					break;
				case C_NW:
					P[i][j] = -(1.0/(2.0*Re))*((U[i][j] + U[i - 1][j] - U[i][j - 1] - U[i - 1][j - 1])/dy + \
						(-V[i][j] - V[i][j - 1] + V[i + 1][j] + V[i + 1][j - 1])/dx);
					break;
				case C_SE:
					P[i][j] = -(1.0/(2.0*Re))*((U[i][j + 1] + U[i - 1][j + 1] - U[i][j] - U[i - 1][j])/dy + \
						(V[i][j] + V[i][j - 1] - V[i - 1][j] - V[i - 1][j - 1])/dx);
					break;
				case C_SW:
					P[i][j] = (1.0/(2.0*Re))*((U[i - 1][j + 1] + U[i][j + 1] - U[i - 1][j] - U[i][j])/dy + \
						(V[i + 1][j] + V[i + 1][j - 1] - V[i][j] - V[i][j - 1])/dx);
					break;
				/* Cells with 2 opposite empty neighbours */
				case C_WE:
					P[i][j] = 0;
					break;
				case C_NS:
					P[i][j] = 0;
					break;
				/* Cells with 3 empty neighbours */
				case C_NWE:
					P[i][j] = 0;
					break;
				case C_SWE:
					P[i][j] = 0;
					break;
				case C_NSW:
					P[i][j] = 0;
					break;
				case C_NSE:
					P[i][j] = 0;
					break;

				/* Isolated cell with 4 empty neighbours */
				case C_NSWE:
					P[i][j] = 0;
					break;
			}
		}
	}
}


inline int is_inside_circle (  	double center_x,
				double center_y,
				double radius,
				double x,
				double y
				)
{
	/* Measure the distance between the point (x,y) and the circle's center and compare it to the radius */
	return ( sqrt( pow((x - center_x), 2) + pow((y - center_y),2) ) <= radius ? 1 : 0 );
}
